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We suggest a simple Gaussian Lagrangian variational scheme for the reduced time-dependent 
quasi-one- and quasi-two-dimensional Gross-Pitaevskii (GP) equations of a dipolar Bose-Einstein 
condensate (BEG) in cigar and disk configurations, respectively. The variational approximation for 
stationary states and breathing oscillation dynamics in reduced dimensions agrees well with the nu- 
merical solution of the GP equation even for moderately large short-range and dipolar nonlinearities. 
The Lagrangian variational scheme also provides much physical insight about soliton formation in 
dipolar BEG. 
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I. INTRODUCTION 

The time-dependent mean-field Gross-Pitaevskii (GP) 
equation can accurately describe many static and dy- 
namic properties of a harmonically trapped Bose- 
Einstein condensate (BEG) [IVH7|. However, the numer- 
ical solution of the three-dimensional (3D) GP equation 
could often be a difficult task due to a large nonlinear 
term [THl IS] . Fortunately, in many experimental situ- 
ations the 3D axially symmetric harmonic trap has ex- 
treme symmetry so that the BEG has either a cigar or a 
disk shape [20 . In these cases the essential statics and 
dynamics of a BEG take place in reduced dimensions. By 
integrating out the unimportant dimensional variable(s), 
reduced GP equations have been derived in lower dimen- 
sions , which give a faithful description of the BEG 
in disk and cigar shapes. For disk and cigar shapes the 
reduced GP equation is written in two (2D) and one di- 
mensions (ID), respectively. The numerical solution of 
such 2D, or ID equation, although simpler than that of 
the original 3D equation, remain complex due to the non- 
linear nature of the GP equation. Hence, for small values 
of the nonlinearity parameter, a Gaussian variational ap- 
proximation is much useful for the solution of these equa- 
tions 

The alkali metal atoms used in early BEG experiments 
have negligible dipole moment. However, most bosonic 
atoms and molecules have large dipole moments and a 
^^Gr USHU], and ^^^T)y BEG with a larger long- 

range dipolar interaction superposed on the short-range 
atomic interaction, has been realized. Other atoms, like 
i66Er [SI 131], and molecules, such as ^Li-^^SCs [37, with 
much larger dipole moment are being considered for BEG 
experiments. A 3D GP equation for a dipolar BEG with 
a nonlocal nonlinear interaction has been suggested [25l 
and successfully used to describe many properties of these 
condensates 

The applicability of the nonlocal GP equation to the 
case of dipolar BEG has been a subject of intensive study 
[25] , After a detailed analysis, You and Yi [36l [S^ con- 
cluded that the GP equation is valid for the dipolar BEG. 



Further support on the validity of this equation came 
from the study of Bortolotti et al. [331 133] . They com- 
pared the solution of the dipolar GP equation with the 
results of diffusive Monte Garlo calculations and found 
good agreement between the two. However, the 3D GP 
equation for a dipolar BEG with the nonlocal dipolar in- 
teraction has a complex structure and its numerical solu- 
tion, involving the Fourier transformation of the dipolar 
nonlinear term to momentum space |35H37| , is even more 
challenging than that of the GP equation of a non-dipolar 
BEG. 

Here we reconsider the dimensional reduction [ISl |35] 
of the GP equation to ID form for cigar-shaped dipolar 
BEG and obtain the precise ID potential with a dipolar 
contact-interaction term. Previous derivations f?SJ |35] 
of the ID reduced equation for dipolar BEG did not in- 
clude the proper contact-interaction term, lacking which 
the ID model will not provide a correct description of 
the full 3D system. We also consider the reduced 2D GP 
equation [171 HE] for a disk-shaped dipolar BEG. Though 
these reduced GP equations for dipolar BEG are compu- 
tationally less expensive than their 3D counterparts, the 
numerical solution procedure remains complicated due to 
repeated forward and backward Fourier transformations 
of the non-local dipolar term. As an alternative, here we 
suggest time-dependent Gaussian Lagrangian variational 
approximation of the ID and 2D reduced equations. A 
direct attempt to derive the variational Lagrangian den- 
sity of the reduced equations is not straightforward due to 
nonlocal integrals with error functions. We present an in- 
direct evaluation of the Lagrangian density avoiding the 
above complex procedure. Thus, the present variational 
approximation involves algebraical quantities without re- 
quiring any Fourier transformation to momentum space. 

In case of dipolar BEG of Gr and Dy atoms we consider 
the numerical solution of the 3D and the reduced ID and 
2D GP equations for cigar and disk shapes to demon- 
strate the appropriateness of the solution of the reduced 
equations. The variational approximation of the reduced 
equations provided results for density, root-mean-square 
(rms) size, chemical potential, and breathing oscillation 
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dynamics in good agreement with the numerical solution 
of the reduced and full 3D GP equations. 



II. ANALYTICAL FORMULATION 

A. 3D GP Equation 

We study a dipolar BEG of N atoms, each of mass m, 
using the dimensionless GP equation 



. ^0(r,t) 

' dt 



1 



+ l^(r) +47ra7V|0(r,i)|^ 



+ N / C/rf<j(r-r')|0(r',t)pdV 



(r,t), (1) 



with dipolar interaction UddCR-) — 3add(l — 3cos^ 6) /R^, 
R = r — r'. Here V{r) is the confining axially symmetric 
harmonic potential, (j){r, t) the wave function at time t 
with normalization J \(j){r,t)\'^dr = 1, a the atomic scat- 
tering length, 9 the angle between R and the polariza- 
tion direction z. The constant add = /xo/2^m/(127r/i^) is 
a length characterizing the strength of dipolar interac- 
tion and its experimental value for ^^Gr is 15ao P5ti28j . 
with oq the Bohr radius, fl the (magnetic) dipole moment 
of a single atom, and fiQ the permeability of free space. 
In equation ([!]) length is measured in units of charac- 
teristic harmonic oscillator length I = y^h/muj, angular 
frequency of trap in units of uj, time t in units of u}~^, 
and energy in units of fku. The axial and radial angular 
frequencies of the trap are flzUj and ^pUj, respectively. 
The dimensionless 3D harmonic trap is 



Vir) = -VL'pP 



(2) 



where r = (p, z) , with p the radial coordinate and z the 
axial coordinate. 

The Lagrangian density of equation ([T]) is given by 

+ 27raiV|0|4+^|0|2y Udd{v-v'mv')\H\'. (3) 
We use the Gaussian ansatz [MJ I55H57] 
(/)(r,i) = ^exp -^r-^ - +iap'^ +il3z^ 



,(4) 

for a variational calculation, where Wp and Wz are time- 
dependent radial and axial widths, and a and /? time- 
dependent phases. The effective Lagrangian L = J Ccfir 
(per particle) becomes 



L = w„a 



2 • I ^zh 



4 2^2 4w2 



N 



+2wla^ + wll3^ + , ^\ , [a - addfi^)] , (5) 



with 



1 + 2«;2 - Sn'^din) 



atanhVl — w„ 
d[K) = ^ — , K = — 

VI - k2 Wz 



(6) 
(7) 



The Euler-Lagrange equations for variational parameters 
Wp,Wz,Oi and /3 yield the following equations for widths 



Wp and Wz 



wp + n Wp = ^ + ^= 5 , 

' wj, V27r w^Wz 



Wz + fllwz = 



J_ ^ 2N [a-addciK)] 
\f2m WpWz 



with 



c(k) = 



2 - - 4k'* + 9K'^d{K) 
(1-k2)2 ' 

1 + 10^2 - 2k^ - 9K^d{K) 



(1 - k;2)2 

The chemical potential p for a stationary state is 



(8) 
(9) 

(10) 
(11) 



P = 



1 , 1 , 2N[a-addf{n)] , , 



p^p ^''z'^z 



2w2 4w2 



27rw2i(;2 



B. ID reduction 



(12) 



For a cigar-shaped dipolar BEG with a strong radial 
trap (fip > rj^) one can write the following effective ID 
equation (details given in Appendix) 



.d4>iD{z,t) 



dt 

°° dk 



92 

y 



2addN 



(f'lDiZjt), 



(13) 



where suj is defined by equation (A9) and dp — l/y/flp 
is the radial harmonic oscillator length. 



To solve equation ( 13 ), we use the Gaussian variational 
ansatz 



-1/4 



exp 



z 

'2u^ 



ifiz^ 



(14) 



From equation ( A2 ) we see that the variational ID ansatz 



( 14 1 corresponds to the following 3D wave function 

^-3/4 



p2 z2 

,w.''''^\ 2dl 2wl 



iPz' 



(15) 



The present variational wave function (15]) is a special 
case of the 3D variational wave function ( 4rwith Wp — dp 
and a = 0. Hence, the ID variational Lagrangian can be 
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written from the 3D Lagrangian (l5|, (using Wp = dp and 



0,) as 



Lid — 



wl(i 1 

H 

2 Awl 

N 



==— [a-a<id/('«o)] ; ^^0 = — 

2TTdlWz Wz 



(16) 



where we have removed the constant terms. This induc- 



tive derivation of the ID Lagrangian ( 16 ) avoids the con- 



struction of Lagrangian density involving error functions 
in the ID potential (AlO) and subsequent integration to 



obtain the Lagrangian. The Euler-Lagrange equation for 
the variational parameter of Lagrangian (161 is 



J_ ^ 2N[a - addciKp)] 



(17) 



The variational chemical potential is given by 



2N[a-aMl{i^o)] , ^ 



Awl 



2-KWzd'^p 



+ 



(18) 



Not only are the above variational results simple and 
yield a good approximation to the ID GP equation, much 
physical insight about the system can be obtained from 
the variational Lagrangian (16). In a quasi-lD system, 



the axial width is much larger than the transverse os- 
cillator length: w^ S> dp. Consequently, — > and 
/(kq) 1- From equation (16), we see that the interac- 



tion term becomes in this limit N{a 
In equation ( 13 1, the dipolar term involves a nonlocal in- 



add)/{V2nd'^pWz 



tegral. However, the variational approximation suggests 
that the effect of the dipolar interaction integral is to 
reduce the contact interaction term in equation ( 13 1 re- 
placing the scattering length a by (a — add)- Immediately, 
one can conclude that the system effectively becomes at- 
tractive for add > a. So one can have the formation 
of bright soliton even for positive (repulsive) scattering 
length a, provided that add > o,- 



C. 2D reduction 

In the disk-shape, with a strong axial trap (fi^ > fip), 
the dipolar BEC is assumed to be in the ground state 

(t){z) = exp(— z^/2d^)/(7rd^)^^^ of the axial trap and the 
wave function 0(r) can be written as [371 EH] 



(TTdl) 



1/4 



cxp 



z 



4>2D{x,y), (19) 



where (f>2D{x,y) is the 2D wave function and d^ = 
■y/l/r^z. Using ansatz (191 in equation ilh, the z de- 
pendence can be integrated out to obtain the following 



effective 2D equation [371 135] 



.d(j)2D{p,t) 

' dt 



AiraddN 



2D\ 



/d^k k d 

exp(ikp • p)n(kp)/i2D(-^) 

"(kp) = y exp(ikp • p) \(j)2D{p)\^dp, 



2'Kdz \/2iTdz 

(f>2Dip,t), 

(20) 
(21) 



where /i2d(C) = 2 — 3y/n£_e^ erfc(^) [l^, kp = {k^^ky) 
and the dipolar term is written in Fourier space. 
To solve equation ([20|, we use the Gaussian ansatz 



<I>2d{p) 



1 



: exp 



2wl 



lap 



(22) 



From equation ( 19 ) we see that the 2D wave function 



( 22 ) corresponds to the following 3D wave function 



(r,t) 



,-3/4 



: exp 



2wl 



z 

2d?, 



lap 



(23) 



The present variational wave function ( 23 ) is a special 
case of the 3D variational wave function ( i^with Wz = dz 
and /3 = 0. Hence, the 2D variational Lagrangian can be 
written from the 3D Lagrangian ^ as 



L 



2D 



Wpa+-^ 



2wl 



2wW 



N 

^w^dz 



(24) 



where we have removed the constant terms. The Euler- 
Lagrange variational equation for width Wp becomes 

J_ ^ N [2a - addgjfi)] 
wl V2^ 



w^dz 



(25) 



The chemical potential p, for a stationary state is 



P- 



2wl 



2N[a - addim 
\/2/Kdzw'l 



(26) 



In a quasi-2D system, the radial width is much larger 
than the axial oscillator length: Wp ^ dz- Conse- 
quently, K — >■ oo and /(k) — >■ —2. From equation (24 1, 



we see that the interaction term becomes in this limit 
N {a + 2add) I {'\/2tt w'^pdz). The variational approximation 
suggests that the effect of the dipolar interaction in equa- 
tion ( 20 ) is to increase the contact interaction term re- 
placing a by (a -|- 2add)- Hence, for positive a, there can- 
not be any bright soliton in 2D, which was found from 
a solution of the 2D GP equation (20) and Bogoliubov 
theory [4^. However, effectively the sign of the dipo- 
lar term in the GP equation can be changed by rotating 
the external field that orients the dipoles much faster 
than any other relevant time scale in the system |50| . In 
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this fashion Nath et al. [SO] suggest changing the dipole 
interaction term by a factor of —1/2, which changes the 
effective scattering length in the Lagrange variational ap- 
proximation to (a — add)i (as discussed in the ID case 
above,) leading to the formation of bright 2D solitons for 
o-dd > o,- These solitons were obtained by Nath et al. 
from a solution of the 2D GP equation (20). 



III. NUMERICAL RESULTS 

We solve the ID, 2D, and 3D GP equations employ- 
ing imaginary- and real-time propagation with Crank- 
Nicolson method TE', TP . The dipolar interaction is eval- 
uated by fast Fourier transform . 

We present results for ^^Cr and ^^''Dy atoms. The ^^Gr 
has a moderate dipole moment with add = 15ao 126H28j . 
while the ^^^Dy atom has a large dipole moment with 
c^dd = 130ao [iniinO]- In both cases we present results 
for dipolar BEG of up to 10,000 atoms for < a < 10 
nm and choose the frequency w such that the oscillator 
length I = 1/im. 

First, we present the results for density profiles ob- 
tained from a solution of the reduced ID and 2D equation 
and compare with the full 3D results. It is known that 
the densities obtained from the reduced equations agree 
well with the full 3D density, as the nonlinearity tends to 
zero and/or the trap asymmetry is extreme [21j . Hence 
in this study we consider a moderately small trap asym- 
metry and a relatively large nonlinearity of experimental 
interest. In the cigar (ID) case we consider ^^Gr atoms 
with a ~ Q nm, and in the disk (2D) case we consider 
i64]3y atoms with a = 6 nm. 
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FIG. 1. Linear density of a cigar- shaped "Cr dipolar BEG of 
1,000 atoms of a = 6 nm, with trap parameters 0.^ = 1 and (a) 
Q,p = 4, and (b) Op = 9 from a numerical (N) solution of the 



3D equation (Tj) and ID equation |13 1, and its variational (V) 
result. Radial density of a disk-shaped ^^''Dy dipolar BEG of 
1,000 atoms of a = 6 nm, with trap parameters Sip = 1 and 
(c) Oz = 4, and (d) Oz = 9 from a nu merical solution of 
the 3D equation ([T| and 2D equation (201, and its variational 
result. 



In Figs. [T](a) and (b), we plot results for linear density 
of a cigar-shaped ^^Gr dipolar BEG of 1,000 atoms as 
calculated from the numerical solution of the 3D equation 
( Tl and the ID equation (13) and its variational result 



(17) for SXz — ^ and i7p = 4 and 9. We find, as the trap 



asymmetry increases by changing f2p from 4 to 9, the 
agreement between 3D and ID models improves. In Figs. 
[T] (c) and (d), we plot results for radial density of a disk- 
shaped ^^^Dy dipolar BEG of 1,000 atoms as calculated 
from the numerical solution of the 3D equation (jlj and 
the 2D equation ( 20 ) and its variational approximation 

We find that, with 
4 to 9, 



(|25j) for f7p = 1 and = 4 and 9. 
the increase of the trap asymmetry from 17^ 
the agreement between the 3D and 2D models enhances. 
In all cases the variational results of the reduced ID and 
2D equations are in good agreement with those of the full 
3D model. 

After having established the appropriateness of the re- 
duced ID and 2D equations in the cigar and disk shapes, 
it is realized that although the numerical solution of these 
reduced GP equations are simpler than that of the full 
3D GP equation, they are still complicated due to the 
presence of the nonlocal dipolar interaction. The vari- 
ational approximation of these equations presented here 
is relatively simple and could be used for approximate 
solution of these equations. Now we test the variational 
results of the reduced ID and 2D equations by comparing 
with the numericalsolution of these equations. 
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FIG. 2. The numerical (N) and variational (V) rms size (p) 
versus scattering length a of a disk-shaped dipolar BEG of 
10,000 (a) ^^Gr and (b) ^'^''Dy atoms for trap parameters Q,p = 
1 and Oz = 4 and 9 from a solution of the reduced 2D GP 
equation (20). The corresponding chemical potential /i in 
these cases for (c) ^^Gr and (d) ^^*Dy atoms. 

In Figs. [2] we present the results for rms size (p) and 
chemical potential /i of a disk-shaped ^■^Gr and ^^''Dy 
dipolar BEG of 10,000 atoms with the trap parameters 
Vlp = 1 and = 4 and 9 for < a < 10 nm as 
calculated from numerical and variational approaches of 
the reduced 2D equation (20). In Figs, [s] we exhibit 
the results for rms size (z) and chemical potential /i of 
a cigar-shaped ^^Cr and ^^'^Dy dipolar BEG of 10,000 
atoms with the trap parameters = 1 and fip = 4 and 
9 for < a < 20 nm as calculated from numerical and 
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variational approaches of the reduced ID equation ( 13 ) 
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FIG. 3. The numerical (N) and variational (V) rms length 
(z) versus scattering length a of a cigar-shaped dipolar BEG 
of 10,000 (a) ^^Gr and (b) ^''^Dy atoms for trap parameters 
= 1 and f2p = 4 and 9 from a solution of the reduced ID 
GP equation ( |13[ ). The corresponding chemical potential fj, in 
these cases for (c) ^'^Gr and (d) ^^*Dy atoms. 

The dipolar interaction changes from strongly attrac- 
tive in the extreme cigar shape (Op fl^) to strongly 
repulsive in the extreme disk shape {flp <C J^z) and its 
effect is minimum (nearly zero) for slightly less than 
ilz ■ In Fig. [2] the dipolar interaction is slightly attrac- 
tive for rtz — 4 and ^Ip — 1. Hence in the absence of 
any short-range interaction (a = 0), the system will col- 
lapse and no stable solution of the GP equation can be 
obtained. For Cr atoms the dipolar interaction is weak, 
and for a > 1 nm, the short-range repulsion for 10,000 
atoms surpluses the dipolar attraction and a stable state 
can be obtained for = 4. For Dy atoms the dipolar 
interaction is stronger, and a stable state can be obtained 
only for a > 2 nm for ft^ — 4. For flz = 9, the dipolar 
interaction for both Cr and Dy atoms are repulsive and 
a stable state is obtained in this case for a > 0. In Fig. [3] 
the dipolar interaction is attractive for both ftp = A and 
9. Hence the dipolar BEC can be stable only for scatter- 
ing length a greater than a critical value. This is why the 
curves in this figure start above this critical value. This 
critical value is larger for Dy atoms and flp — 9 com- 
pared to that of Cr atoms and fip = 4 as can be found 
in Fig. [3j As there is no real collapse in ID models with 
cubic nonlinearity; for confirming the collapse correctly 
one must solve the full 3D GP equation. 

Next we study, by numerical and variational solutions 
of the reduced ID and 2D equations, the dynamics of 
breathing oscillation of the four dipolar BEC of cigar- and 
disk-shaped Cr and Dy atoms shown in Figs. [T] started 
by a small change of the scattering length. This can 
be implemented experimentally by a Feshbach resonance 
[5Tj. In Fig. [4] this dynamics is shown for a cigar-shaped 
Cr dipolar BEC of 1,000 atoms for (a) r2p = 4 and (b) 
rip = 9 from a solution of the reduced ID equation, and 
for a disk-shaped Dy dipolar BEC of 1,000 atoms for (c) 
= 4 and (d) fi^ = 9 from a solution of the reduced 2D 
equation. In these figures we also show the results from 



FIG. 4. The rms sizes (z) of a cigar-shaped Gr dipolar BEG 
of 1,000 atoms versus time t for (a) ilp = 4 and (b) 9 from a 
numerical (N) and variational (V) results of the reduced ID 
equation. The rms sizes (p) of a disk-shaped Dy dipolar BEG 
of 1,000 atoms versus time t for (c) il^ = 4 and (d) 9. The 
oscillation was initiated by jumping the scattering length a 
from 6 nm to 6.15 nm for Gr (6.3 nm for Dy) at t = from a 
solution of the reduced 2D equation. 



a numerical solution of the 3D Eq. ([T]). The agreement 
between the numerical and variational results is good in 
all cases. We also calculated the angular frequencies of 
these oscillations. In case of Cr in Figs. |4] (a) and (b), the 
axial frequencies are 1.75 (variational, ID), 1.76 (numer- 
ical, ID) and 1.63 (numerical, 3D), and in case of Dy in 
Figs. [4](c) and (d), the radial frequencies are 1.93 (vari- 
ational, 2D), 1.89 (numerical, 2D) and 1.76 (numerical, 
3D). For quasi- linear systems, these angular frequencies 
are expected to be 2 [52]. The deviation from this value 
is due to the large nonlinearity of the dipolar BECs con- 
sidered here. 



IV. CONCLUSION 

The usual GP equation provides a good description of 
statics and dynamics of a normal nondipolar BEC. For 
a dipolar BEC the numerical solution of the GP equa- 
tion is a difficult task due to the nonlocal dipolar inter- 
action. For a cigar- and disk-shaped dipolar BEC, the 
reduced ID and 2D equations provide an alternative to 
the full 3D equation. Nevertheless, the solution of these 
reduced equations is also challenging involving Fourier 
and inverse Fourier transformations. As an alternative, 
we suggest a time-dependent variational scheme for these 
reduced equations, not requiring any Fourier transfor- 
mation. The variational approximation of these reduced 
equations provides results for stationary cigar- and disk- 
shaped dipolar BEC as well as for breathing oscillation of 
the same in good agreement with the numerical solution 
of the respective GP equations. This is illustrated for 
large Cr and Dy dipolar BECs of 10,000 atoms and large 
atomic scattering lengths a up to 20 nm. We also study 
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the breathing oscillation of a bright soliton of 1,000 Cr 
atoms using the numerical solution of the 3D equation 
as well as the numerical and variational approaches to 
the ID equation. A typical dipolar BEC considered here 
corresponds to a large short-range cubic nonlinearity of 
about inaN « 1250 for a = 10 nm and N = 10, 000 and 
a large dipolar nonlinearity of '^Traj^^ w 865 for Dy atoms 
for add — ISOflQ and N = 10,000. The variational ap- 
proximations considered here provided good results for 
such large nonlinearities and should be useful for analyz- 
ing the statics and dynamics of realistic dipolar BECs 
under appropriate experimental conditions. 
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— oo ^ —oo 
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where the ID potential in Fourier space is 
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Jo 
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Appendix A: ID reduction 



For a cigar-shaped dipolar BEC with a strong radial 
trap (r2p > ^z)^ we assume that in the radial direction 
the dipolar BEC is confined in the ground state 



(Al) 



of the transverse trap and the wave function ^(r) 
0i£)(z) X(f>{p) can be written as [iSlli^ 



(r) 



: exp 



4>iD{zy, npdi - 1, (A2) 



where dp is the radial harmonic oscillator length. 
The contribution of the dipole potential to energy is 



N 
~2 
1 N 



Hdd^j I d-'r I d\'n{T)Udd{v-v')n{T'), (A3) 



2 (27r)3 



d3fcn(k);7dd(k)n(-k), 



(A4) 



where ri(r) = |(/'(r)|2 is the density and in Eq. (A4) we 
used a convolution of the respective variables to Fourier 
space and where tilde denotes Fourier transformations: 



sid(C) = 



/"OO 

/ du 
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(A9) 



The ID potential in configuration space is 
1 1"°° 

U'dF(Z) = / dkze^'^'^VMkz) 



{V2dpf 



-d{Vi) + 2Vt- V^(l + 2^)e*erfc(Vt) 



(AlO) 



Similar, but not 



where t = [Z/{V2dp)]^, Z = \z - 
identical, ID reduced potential was derived in [151 
where the (5- function term was absent. To derive the 
effective ID equation for the cigar-shaped dipolar BEC, 
we substitute the ansatz ( |A2| ) in Eq. ([ij , multiply by the 
ground-state wave function (/'(p) and integrate in p to get 
the ID equation 
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